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We examine the time-dependent dynamics of ion crystals in radiofrequency traps. 
The problem of stable trapping of general three-dimensional crystals is considered 
and the validity of the pseudopotential approximation is discussed. We derive analyt- 
ically the micromotion amplitude of the ions, rigorously proving well-known exper- 
imental observations. We use a method of infinite determinants to find the modes 
which diagonalize the linearized time-dependent dynamical problem. This allows 
obtaining explicitly the ('Floquet-Lyapunov') transformation to coordinates of de- 
coupled linear oscillators. We demonstrate the utility of the method by analyzing 
the modes of a small 'peculiar' crystal in a linear Paul trap. The calculations can be 
readily generalized to multispecies ion crystals in general multipole traps, and time- 
dependent quantum wavefunctions of ion oscillations in such traps can be obtained. 
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I. INTRODUCTION 

The trapping of charged particles by electromagnetic fields is an essential tool for many 
investigations within physics and chemistry as well as for studies of bio- molecules 

pi, [a]. In the Paul trap, charged particles are confined by oscillating multipole radiofrequency 
(rf) fields jij]. Different types of Paul traps have been proposed and constructed over the 
years, and an increasing number of experimental and theoretical work is dedicated to the 
improvement of these devices (see, e.g. Q). With the advent of laser-cooling techniques 



a lot of effort is focused on the study of crystalline forms of trapped ions [9M37] . Single 



trapped ions have been cooled to the quantum ground state of motion [38|, |39| , and crystals of 



a few ions have been cooled to the ground state in at least a few of the motional modes 40 



Long-range order has been observed with large structures of thousands of ions in Penning 



2i|, |22j and Paul traps [37], |41 



In this contribution we analyze the dynamics of crystals of charged particles confined in 
rf traps. In Sec. [TT] we review some previous results on the trapping of ions in rf traps. In 
Sec. lIIII we discuss the linear stability of ion crystals in general, and consider limits of validity 
of the pseudopotential (time-independent) approximation, in relation to symmetries of the 
trapped crystal. We derive the micromotion amplitude of ions in a general periodic solution 
in a Paul trap, and relate our findings to recent experimental results. We then expand the 
motion of the trapped ions in coordinates of small displacements about the periodic solution, 
which is the dynamic equivalent of a minimum of the trapping potential. This expansion is 
completely general and can be applied to any periodic trapping field. In Sec. [TV] we solve 
the coupled motion of the ions in the time- dependent potential. We find the modes which 
diagonalize the dynamical problem and obtain explicitly a time-dependent transformation to 
coordinates in which the motion is that of decoupled linear oscillators. Using this expansion 
the exact time- dependent wavefunctions of ions in rf traps can be obtained. We present a 
numerical study of a small crystal in a linear Paul trap in Sec. |Vl and conclude in Sec. IVII 
with comments on further possible applications and research. 



II. OVERVIEW OF PAUL TRAPPING 



The first crystallization of a cloud of charged aluminum microparticles has been reported 
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Wigner crystals of 2 to approximately 100 trapped ions in a Paul trap were re- 
ported in [9|, llO] and further investigated in 43[ both experimentally and numerically. The 
simulations included the rf trapping, Coulomb interaction, laser cooling and random noise. 
Depending on trap parameters, the ions have been found to equilibrate either as an appar- 
ently chaotic cloud or in an ordered structure. The latter is defined as the 'crystal' solution 
when it is a simple limit cycle, with the ions oscillating at the rf frequency about well-defined 
average points. The transition between the two phases has been investigated, it was shown 



that both phases can coexist and hysteresis in the transition has been observed 43 



The motion of two ions in the Paul trap has been investigated in detail in various publica- 
In addition to the aforementioned phases, frequency-locked periodic attractors 



tions 44-53 



(where the nonlinearity pulls the motional frequencies into integral fractions of the external 
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rf frequency) were found in numerical simulations and experiments. These solutions are dif- 
ferent from the crystal in that the ions move in extended (closed) orbits in the trap, whose 
period is a given multiple of the rf period. However many of these frequency locked solutions 
are unstable, especially those of a large period, and perturbations such as those coming from 
the nonlinearity of the laser-cooling mechanism, tend to destroy them. 

Despite the large amplitude motion, the frequency-locked solutions, being periodic, are 
of course not chaotic. However, even at the presence of cooling, some solutions in the 
ion trap may behave chaotically for exponentially long times. The authors of 54] suggest 
that, eventually, all trajectories settle at frequency locked attractors, at least for two ions 
at a = 0. Numerical simulations and experiments with more ions suggest that in general, 
different types of solutions - of chaotic and of long-range order nature, can coexist at the 



same parameter values 43 



A further important property of two-ion crystals was discovered in [55(. The entire first 
Mathieu stability zone (which corresponds to stable trapping of one ion), was mapped using 
numerical simulations to determine the stability of the two-ion crystal. It turned out that the 
Coulomb interaction destabilizes completely the two-ion crystal at some parameter areas. 



About a year later, it was independently discovered in [56|, |57| and [58], that two ions in a 
hyperbolic Paul trap may crystallize in a 'peculiar' crystal. 

The reason that such crystals were termed 'peculiar', is that for a two-ion crystal with no 
angular momentum about the axial x-axis, there are only two possibilities in the harmonic 
pseudopotential approximation (excluding the case of degenerate secular frequencies). The 
crystal can either align with the x-axis or lie in the y-z plane. In the latter case we may 
assume that the crystal is aligned with the |/-axis, and in both cases the z coordinate can 
be eliminated. Therefore a crystal which is at angle to both reduced axes x and y, i.e. one 
which forms an angle with both the axial axis and the plane of symmetry, is peculiar. 



In 5a, |57|], an improved pseudopotential approximation was derived for two ions in the 



hyperbolic Paul trap (and in [59] for the linear Paul trap). This nonlinear pseudopotential 
can be used to derive approximately the areas of stability in parameter space of the axial 
and radial two-ion crystal, and the orientation of the peculiar crystal. However, even the 
improved nonlinear pseudopotential cannot reproduce the areas of no stable crystal. 

As mentioned above, the crystal solution in all of these works, is taken to be a periodic 
solution of the equations of motion with a period equal to that of the rf potential. This is a 
limit cycle of the equations, the simplest of the frequency locked solutions mentioned above, 
and it is an attractive one when cooling is present and it is stable. Its stability is analyzed 



m 



58[ by looking at the Poincare map associated with this solution. The Poincare map is 



a mapping of the phase space into itself, in which every inital point {x (0) ,p (0)} is evolved 
according to the equations of motion and mapped to {x (T) ,p(T)} after one rf period T. 
The crystal is by definition a fixed point of the Poincare map, and its linear stability is 
determined by the linearized dynamics in its vicinity, specifically by the eigenvalues of the 
so-called monodromy matrix. The linearization takes into account the leading (harmonic) 
coupling between the ions, expanded about the periodic solution. Further analysis (both 
linear and nonlinear) of the periodic orbit of two ions in a Paul trap appears in 
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III. PAUL TRAPS 

A. The Trapping Potential 

We start with the potential energy of identical ions in the most general single, nonseg- 
mented quadrupole rf trap. By choosing a specific frequency Co, which is characteristic of 
the secular frequencies in the problem, and measuring time and distances in units of 1/ui 

1 /3 

and d = (e 2 /muj 2 ) , respectively (with m and e the ion's mass and charge), we write the 
nondimensional potential as 



N 1 1 II -i 

V = Vtrap + Coulomb = }2 9 + ^ + + Yl 2 11^ ~ ^ ' ^ 



where Ri = {xi,yi,Zi} is the vector coordinate of ion i, the trapping terms are given by 
A a = -j- (a Q — 2q a cos fit), a E {x, y, z}, with a a and q a being the nondimensional Mathieu 



parameters of the respective coordinates [64|, and being the rf frequency (in units of u). 
Regarding the trapping parameters, the Laplace equation implies that 

5> Q = 5>a = 0. (2) 



Two commonly used types of traps are the hyperbolic [65j and the linear [64] Paul traps. 
Taking the the axial direction, the hyperbolic trap can be described by setting 

a y = a z = a, a x = —2a and q y = q z = q, q x = —2q, while in the linear trap we have 
a y = a z = a, a x = —2a and q y = —q z = q, q x = (so a must be negative to obtain 
stable trapping). In the next subsection we will relate to symmetries of these two trapping 
geometries, but for now we keep the discussion completely general. 

The trapping potential Vt rfl . v w hen considered alone, gives rise to decoupled Mathieu equa- 
tions in each ion coordinate [66J. The corresponding characteristic exponents are /3 a (a a , q a ) 
and the secular frequencies are then u a = [3 a ^. In these units, V tra , p and Vcouiomb are both 
of order unity for a crystal, when the trapping balances the Coulomb repulsion. These units 
allow naturally the introduction of the parameter 

e = 4/fi 2 . (3) 

If u is the (dimensional) secular frequency along some specific axis a, we have u a = 1 and 



e = f3 2 - Using the familiar lowest order approximation 64j. 8„ \J a a + q%j2, we see that 
the limit e — > is equivalent a — > 0, q 2 — > 0. This parameter will be used in the following. 

When the Mathieu parameters belong to a single-particle stability-zone, the potential of 
eq. (00) allows for bounded motion and large ion clouds can be trapped for extraordinary 
long times without cooling. This is given by the limit Ri — > oo in eq. ([I]), which means that 
ions are not crystallized, rather their motion is that of decoupled trapped particles. Despite 
the large amplitude motion, at low enough density the ions hardly interact. 

However, stable single-particle parameters do not guarantee that a stable crystal solution 
exists, even with two ions, as discussed extensively in Sec. [TT1 which has also been recently 
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investigated with crystals consisting of upto several hundreds of ions 67], |68(. We recall 



that a stable crystal solution is considered as a periodic solution with the same period of 
the rf, which is linearly stable under small perturbations. The existence of a stable crystal 
solution is a property of the fully nonlinear and time- dependent problem. There may also 
exist stable multi-periodic solutions (with a period which is some large, arbitrary multiple 
of the rf period), the frequency- locked solutions described in Sec. flU 

In experiments, even at the presence of cooling, crystals may exist only metastably, 
changing completely after a given time or when parts of the crystal are moving with respect 



to the rest 37]. In this case, if the timescale of the change of the crystal is long enough, 
the analysis in the following subsections, using modes expanded about the (quasi-) periodic 
solution may still be useful, even though at least one mode would be unstable. 



B. Linearization using the Pseudopotential Modes 



In this subsection we linearize the crystal motion using the pseudopotential modes. There- 
fore this treatment assumes that a crystal solution exists, and is close to the crystal of the 
pseudopotential. As known from the example of two-ions in a hyperbolic Paul trap, this 
assumption may not hold, even for arbitrarily small values of a and q (where the crystal is 
'peculiar' 56l-l58|). The linearized secular modes of a one-dimensional Coulomb chain were 



treated in detail in 169-71 



Let us define the squared ratio of secular radial and axial frequencies 7 y = Uy/u^, 7 Z = 
u)z/oJx and j x = ojI/^I = 1> thus choosing u as the axial trapping frequency, so we can 
rewrite the potential as the sum V = Vp Seu d + Vcouiomb + Kfj 

N * „ N 



i,a i^tj i,a 



(4) 



We expand about the minimum-configuration locations, j> obtained from the secular 
part of V in eq. fll]), Vp Seu do + ^Coulomb, and change to the normal modes 0j by setting 

3iV 



R ha (t) = Rl a + Dl«®j W 



(5) 



where D J ia is the matrix of normal mode vectors, with rows indexed by the iV ion numbers i 
and 3 directions a, and columns by the 3iV normal mode numbers j. Writing the potential 
in terms of the normal modes, V = Vharmonic + V T { + Vnoniinear, we keep only the first two 
terms to get 

2 



N /37V 

^=£^ 2 +£^(a«-7 Q ) Ui+y: 

i i,a \ j 



+ 



and the linearized equation of motion (e.o.m) derived from eq. ([6]) is 

dV r{ 



<96 r 



N / 37V \ 

- £ (A Q - la ) D™ a Rl a + ^ D^Qj 

i,a \ j J 



(6) 



(7) 
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Before treating the coupled system described by eq. ((Tj), we first note that when the rf 
trapping potential is symmetric with respect to the axes of symmetry of the crystal, the 
equations of motion are in fact diagonal. For the following few paragraphs, let us therefore 
limit the discussion to the hyperbolic and the linear Paul traps, whose parameters were 
described following eq. (j2J). Given these two trapping geometries, the current expansion will 
be diagonal if the crystal is two-dimensional and lies in the plane of cylindrical symmtery of 
a hyperbolic trap, or the crystal is one- dimensional and aligned with one of the trap axes, in 
either the hyperbolic or the linear trap. In that case, we have R® a = for the coordinates 
transverse to the crystal, and also the normal modes decouple in these directions from the 
modes tangential to the crystal (i.e. D{ a is divided into blocks in the index a). 

We then can use two relations which hold in the plane of the crystal or along its axis, 

/] D i?a D i& = *W) /] D i!kRla = CbSrnb, (8) 



where hereafter, a runs on the symmetry directions , £& = y Yli & an d b denotes 

the breathing mode. The identity on the left in eq. ([8]) is the completeness of the normal 
modes along the symmetry directions. To get the second identity we use the fact that 
in a harmonic trap, the breathing mode vector is exactly proportional to the minimum- 



configuration locations R® a of the pseudopotential (see, e.g. App. B of [l5|), so that 
D\~ = R® & /£b, and therefore the vector i?° a is orthogonal to all other normal mode vectors. 
By using eq. (jSJ) we can replace the r.h.s in eq. (j7|) with the simple expression 

dV 

- ~a7T~~~ = ~ (Afi - la) (ZbSmb + e m ) • (9) 

We find that the equations of motion for modes in the crystal plane or along its axis are, 
after multiplying by e = 4/Q 2 of eq. fl3]) and rescaling t — > Qt/2, 

B m + [e (co 2 m - 7) + (a - 2q cos 2t)] 6 m = [ery - (a - 2q cos 2t)\ £ b 5 mb (10) 

where 7, a and q are defined along the symmetry axes. Similar equations hold in the direc- 
tions transverse to the crystal, without the inhomogenous r.h.s and with the corresponding 
7, a and q. 

Eq. fllOp shows explicitly how the isotropy of the potential along the axes of symmtery 
of the crystal allows to decouple the equations of the modes. This decoupling puts the 
conditions for linear stability of the crystal in terms of diagonal Mathieu equations for each 
of the modes. In addition, with this symmetry the only mode with an inhomogenous r.h.s 
is the breathing mode, i.e. it is the only mode on which the rf potential acts as a driving 
force. This driven motion is in fact the 7r-periodic motion of the crystal as whole, about the 
static minimum-configuration locations of the pseudopotential, j-Rfj. 

Returning to the case of a general trap, in order to put eq. (JTj) in a simple form we 
multiply it by e = 4/f2 2 of eq. ([3]) and rescale t — y Qt/2. Defining the two matrices 

A mj = eU 2 J m3 + J2 D Z D U ( a * ~ , Qmj = D Z D U^ (11) 
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and the vectors 

G m = - A> ( a a - e7o) Rl a , F m = ^2 D™ a q a R° i al (12) 

i,a i,a 

we rewrite eq. © in vector notation as 

6 + [A-2Qcos2t]Q = G + 2Fcos2t. (13) 

We solve the homogeneous l.h.s of eq. ([TBI in Sec. [TV] Since the pseudopotential modes 
are expanded aboout static configuration points, eq. (fTB"]) is an inhomogenous equation with 
a 7r-periodic r.h.s, which has a unique 7r-periodic solution (except for possibly a region of 
measure-zero in a a , q a space). This periodic solution of driven motion of the normal modes, 
corresponds to the exact ^-periodic solution which defines the crystal in the rf trapping 



potential. Details of the solution of the inhomogeneous equation can be found in [72 . 

Eq. (TIB"]) shows that in general, the rf couples the pseudopotential normal modes and also 
acts as a driving force. Under this coupling, the true modes of oscillation of the system 
may in general have different frequencies (and even lose stability), and different oscillation 
directions than the pseudopotential modes upon which this expansion is based. Indeed, the 
linearization starting from the pseudopotential approximation may not be adequate in the 
general case. We investigate this point further in the following two subsections. 



C. The Periodic Crystal Solution 

We now abandon the pseudopotential approximation and turn to studying the time- 
dependent potential directly. In this subsection we derive analytically the micromotion 
amplitude of the ions in typical crystals in Paul traps. The e.o.m for the ion coordinates, 
derived from eq. (CEJ) after rescaling by t — > Qt/2, is 



N -3 

Ri,a + (o« _ 2g a cos 2t) R i: a - e ^ LRj - Rj ^ 



3=1 



R% 



R 



0. 



(14) 



Eq. (114j) has ^-periodic coefficients and is time-reversal invariant. We assume the existence 
of a crystal in the sense of Sec. UH i.e. a 7r-periodic and time-reversal invariant solution, 
which obtains the genreal form 



(15) 



In this form, the average ion location is B i (which is different from RJ (t = 0)). 

We wish now to see what can be said about RJ in typical Paul trapping experiments. 
We first define the dynamic matrix 



Ri{t) 



5ij ^ R{ Rn 



m 



-(1-5. 



13) 



R,; — Rn 



(16) 
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and write using eq. (TT5T) its Taylor and Fourier expansion around j-Bo,i — as 

C« = Co,ii ({^o,i}) + G 2>ij ({S 0ii } , {s 2ii }) (e 24 * + e" 24 *) + ... (17) 



Substituting the solution ansatz eq. (fT5|) into eq. (|T4|) we get 



(a Q — (2n) 2 ) B2n,i,a — Qa (^B2n-2,i,a + B2n+2,i,a^j ~ C GijB 2n j : , 



e l2nt = 0. (18) 



We replace Gjj ^|-Ri(t)|j in eq. ( IT8|) by its leading order, time-independent term from 
eq. ( fTTl) . Go,ii = Gjj ^ji?o,jjj, and require that the above relation holds for every t. We 

get, by using B 2n = B-2n and neglecting B in « (which is of the next order in g a ), the 
equation coming from the n = term, 

— 2q , Q -B 2j i ja — e Gq^Bqj^ = 0, (19) 
i 

and the coupled equation from the coefficient of e 2lt , 

(a a — 4) B 2t i t a — q a Bo,i,a — e ^ G 0ti jB 2 ,j,a = 0. (20) 

j 

This sytem of equations can be seen as a linear homogenous system in the 2N variables 
|-B 0) j )a , B 2> i >a j, i = 1,...,N, for fixed a, since the equations are diagonal along the three 
different axes. Of course, this is not really a linear system since the matrix Gq^ depends on 
B 0) i, but since we are not trying to actually solve the system, this will not matter. Let us 
fix the index a to some axis (supressing it in the following), and define the iV-component 
vectors u 0) i = B 0)i>a and u 2;i = B 2>i>OL . Then the above system can be written in block- matrix 
form 

a — eG —2q \ / u 
-q -4 )\u 2 



0, (21) 



where we have neglected a — eGo as compared with —4 in the lower- right block of eq. (T2TT) . 
Since we assume that the system has a solution (which approximates the 7r-periodic crystal), 
the above matrix must be singular. We can expand its determinant, 

= det [a - eG ~ (~2q) (-4) -1 (-<?)] = det [a - eG + g 2 /2] . (22) 

Taking u to be the vector from the kernel of (a — eGo + <? 2 /2) which obeys (a — eGo) = 
— (q 2 /2) uo, we find that the solution of eq. ( 12~TT) is u 2 = — (g/4) u Q . 

We therefore have obtained that in a general quadrupole trap, in the 7r-periodic crystal 
solution of eq. (|T5|) . every ion coordinate obeys (at least to leading order in a Q /4, g a /4, e a /4), 

B 2i i >a ~ — yBo,i,a, (23) 
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i.e. that the micromotion amplitude in each coordinate is q a /2 of the respective average 
position, and at 7r phase with respect to the rf drive. In the hyperbolic trap the corresponding 
motion has been imaged as early as in {4^. In simulations of a generic trap (with different 
q and a parameters for the three axes), eq. (1231 seems to hold accurately to within a few 
percent, for q up to ~ 0.7, which is consistent with a deviation of order ((fo/4) . 

The relation in eq. ( 123]) loses its accuracy when either -E>o,i,a <C 1 or q a <C 1. In the former 
case, for an ion near the zero of one of the rf axes, the corresponding micromotion amplitude 
in fact seems in the cases we have checked, to be lower (this is similar to a single trapped 
ion at the origin of a Paul trap, for which the unique 7r-periodic solution is the trivial one). 

For the linear Paul trap case of q a <C 1 along the axial axis, the first-order expression in 
eq. (123]) loses its meaning. We must therefore add to eq. (l20~l) the second term in eq. ( fl7j) . 
G 2 ,ij, which rotates at the frequency e 2tt . Setting q = 0, eq. (121]) is replaced with 



a — eG \ / u Q 
-eG 2 -4 J I u 2 



0. 



(24) 



We now examine what can be deduced about G 2 . The off-diagonal elements are in general 



G 2} ij — 3 



Bo,i — -So, j J ^ (j3o,i,a — Bqj^ \B2,i,a — B 2! j,aj , % 7^ j. (25) 



In particular, using eq. ( 123]) with q y 
trap 



-q z = q and q x = 0, we get that for the linear Paul 



G 



2,ij 



Bqa — B t 



Bo,i,y ~ Boj^y 



Q.i.z 



B, 



0,j,z 



off. 



i + j- 
(26) 



Since the diagonal elements are the negative of the row sum, as in eq. ([16]) . G 



2,ii — 

~ Ylm^i ^2,im, we immediately find that for a crystal which is invariant (up to a permu- 
tation of the ions) under y -h- ±z, G 2) u = + O (g 2 /4 2 ). In fact, for the linear Paul trap, 
the e.o.m, eq. (114"]) is invariant under y <H- —z,t — > t + 7r/2, so given one crystal solution, 
there is also a solution related by this transformation. Depending on the number of ions 
and trapping parameters, both solutions may actually be the same crystal. As the crystal 
size grows, by the same symmetry arguments, G 2jii — > + O (g 2 /4 2 ). 

The second row of eq. (1241) gives B 2 a,x = ~ e /^J2^^,ijBo,j,x, and we argue that by the 
above symmtery arguments, for a typical symmetric or large crystal in a linear Paul trap, 
when the crystal configuration is also symmetric under x -H- —x, the first order terms in the 
above summation will cancel, and (using e ~ q/2) 



B 



2,i,x 



I \- fa ) Do.. 



(27) 



If the symmetry alluded to above does not hold, eq. (J27j) should be replaced with an expres- 
sion which is one order less, namely B 2 ^ x — O (| • |J B Q i x . 

Eq. ( 127]) explains why in the linear Paul trap, there is essentially no micromotion exci- 
tation along the axial direction despite the strong Coulomb interaction. In addition, since 
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q y = —q z in the linear trap, the oscillation described by eq. (|23|) is exactly the (2,2) mode of 
cold-fluid theory [14| which has been observed both in simulations 20 j and recently discussed 



in connection with experiments 29) . In fact, in simulations we have performed (Sec. |V|) . 
eq. (1231) seems to hold extremely accurately radially (to within half a percent), and eq. (ET 
gives an accurate estimate for the axial micromotion amplitude (also matching e.g., (20 



D. Linearization about the Periodic Crystal, and the Pseudopotential Limit 

In this subsection we expand the modes of oscillation of the ions about the periodic crystal 
solution, eq. ( TT5l) . The following expansion is applicable to general crystal configurations, 
does not rely on the pseudopotential modes or on a small parameter, and can be readily 
generalized to settings not considered in this paper, e.g. segmented traps and higher-order 
multipole traps. We will discuss its relation to the expansion of Sec. IIII B\ which used the 
pseudopotential modes. 

Expanding the potential of eq. ([1]) about the time- dependent solution |_R^ (t)|, in the 

coordinates of small oscillations r i<a = R ia — RJ a (t) , and defining the 7r-periodic matrix 



f) 2 V 

K°nt) ( " ui * ,ui1 ' 



" "' dR it(T dR h7 



(28) 



we find the linearized e.o.m 

TV + (a* - 2q a cos 2t) r i>a + e ^ K% (t) r hT = 0. (29) 



This equation is a linearly coupled system for the ion coordinates with 7r-periodic coeffi- 
cients. We see that the interaction term in eq. (I2"9"j) is of the same order (in e, or equivalently 
a and q 2 ) as the diagonal term. Expanding the matrix K^J it) in a Fourier series in the form 

K°l = (K )°J -2(K 2 )%cos2t-..., (30) 

we first define the two matrices 

A1J = hA^c + e (KoYJ , Q1J = 5 l3 5 aT q a + e (K 2 )% . (31) 

We can now switch notation to a simpler one with dynamical variables u m whose single 
index m stands for both indices {i, a} of rj jCr , with the corresponding indices replaced in 
eq. fl3Tl) . and rewrite eq. fT29|) in vector notation, with only the two leading harmonics of the 
Fourier series, as 

[A-2Qcos2t}u = 0. (32) 

This equation describes linealrized coupled perturbations about the 7r-periodic solution 
which defines the crystal in the time-dependent potential. In Sec. IIVI we will solve this 
coupled system and find its decoupled modes of oscillation. 

Eq. fl32|) is seen to be identical in form to the l.h.s of eq. fflBl . which is given in terms of the 
pseudopotential modes. However, the current expansion is based on the exact force acting 
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between the ions during their motion along the periodic trajectory of the full potential. 
Higher harmonics from the Fourier series of this force can be added, and in Sec. |V]we will 
in fact include the next harmonic (cos At), to obtain very accurate results for the modes. 
When the nonlinear motion is not exactly 7r-periodic, some of the linearized modes of the 
above expansion may not be oscillatory (as detailed in Sec. IIV Aj) . which describes (locally) 
the aperiodic motion of the crystal. If this change of the crystal is slow, the linearization of 
eq. ( l29i) can still be useful. 

Considering the limit of taking the Mathieu parameters to zero, this does not guarantee 
that the crystal will be (linearly) stable. If there are unstable modes, they may remain 
unstable even as a, q 2 , e — > 0. As for approaching the pseudopotential crystal in this limit, 
if the average ion locations j-Bo,ij of the periodic solution tend to the pseudopotential 
minimum locations, then the limit of the pseudopotential modes is regained on the l.h.s 
of eq. ( 1291) . Characterizing the conditions for this to occur requires further investigations, 
however for simple cases the results of Sec. IIII CI may provide some hints. 

In particular, if the rf potential is isotropic with respect to the configuration (as discussed 
in Sec. IIII Bl for specific cases), and we neglect K 2 in eq. ( 13T|) . the matrix Q becomes pro- 
portional to the identity. Then we can diagonalize eq. fl32|) by an orthogonal transformation 
(which is exactly the transformation to the normal modes) and obtain eq. (fTOj) . There will 
now not be an inhomogeneous r.h.s, since the periodic crystal motion is already accounted 
for. At the presence of driven breathing oscillations, the ions feel stronger repulsive nonlin- 
earity when they oscillate radially inwards, so we may expect the real crystal to be a bit more 
expanded (with the ion distances somewhat larger) as compared with the pseudopotential 
approximation. 

We note that for an axial chain of ions in the linear Paul trap, indeed K 2 ~ even if 
there is a small rf leak in axial direction, or when the ions do not lie exactly along the rf null 
axis, and eq. ( ITUl) is almost exact. For the axial modes, the correction to the pseudopotential 
modes will be small. For the radial modes, the correction will typically be small, provided 
that the rf frequency is large compared with the radial frequencies. However, for calculation 
of small effects which might be of importance in high-accuracy experiments such as for 



quantum information processing [64J, the full time-dependent solution must be used. 

As a final note, a generalization of the isotropy of the rf potential with respect to the 
crystal configuration, would be the existence of a constant orthogonal transformation which 
diagonalizes eq. (132]) . and thus decouples the modes in the non-isotropic case. Such a trans- 
formation would simultaneously diagonalize the matrices A and Q (which must commute), 
and does not exist in the general case. 



IV. SOLUTION OF THE LINEAR EQUATIONS 
A. The Floquet Problem 

Eq. ( 132]) is a linear differential equation with periodic coefficients and therefore amenable 



to treatment using Floquet theory, which we briefly review here. For more details see 72 



and references within. For the Newtonian problem with / degrees of freedom (/ = "SN 
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for iV ions in three-dimensions), the corresponding Floquet problem is stated in terms of 
coordinates in 2/-dimensional phase space by the definitions 

*=(i). ^^(-(^-"qcos^o)' (33) 

where 1/ is the /-dimensional identity matrix. The e.o.m is written in standard form as 

<p = n(t)<p. (34) 

In the following, an /-dimensional vector u will be denoted by a lower case Latin letter 
with an arrow, /-dimensional matrices will be denoted by capital Latin letters iff). A 
2/-dimensional vector <p will be denoted by a lower case Greek letter. Capital Greek letters 
(unitalicized) will denote 2/-dimensional matrices (II, B). 

A fundamental matrix solution to eq. f )34|) has 2/ linearly independent column solutions 
and obeys the matrix equation $ (t) = II (t) $(;£). A fundamental matrix solution which 
equals the identity matrix at t — 0, i.e. obeys $ (0) = I2/, is known as the matrizant 
of eq. ( )34l) (and is obviously unique). The matrizant can always be written in the form 
$ (t) = r (t) e^T' 1 (0) where T(t + T) = T (t) is periodic with the period T of II (t), and 
the constant matrix B is diagonal, with entries known as the characteristic exponents of the 
Floquet problem, 

B = diag {iPi, i/3 2 /} • (35) 

The time-dependent linear coordinate change, known as the 'Floquet-Lyapunov' transfor- 
mation, defined by 

0(t)=r(t) X (t), (36) 
transforms the e.o.m eq. (|34|) into the time-independent diagonal form 

X = B X , (37) 

whose solutions are the Floquet modes, 

xAt) = Xu(0)e^ t . (38) 

We note that in the general case, T mixes the coordinates u and their derivatives, and 
in addition is not unitary, although certain highly symmetric cases, e.g., as discussed in 
Sees. IIIIBI and IIII D\ when the Mathieu equations decouple completely, are an exception. 



B. Solution Using an Expansion in Infinite Continued Matrix Inversions 



We expand the solutions of the e.o.m (l32l) b y u sing an analytical expansion which utilizes 



infinite-continued matrix inversions 74 



sec 



72 for details and further references. This 



expansion allows to obtain the frequencies and the coefficients of the solution vectors, to ar- 
bitrary precision. The solution is given in a form which is immediately suitable for obtaining 
the Floquet-Lyapunov transformation, as will be shown in the next subsection. 
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We seek for solutions of eq. ( 132]) . in the form of a sum of two linearly independent complex 
columns vectors, 



u 



[be l(2n+ ^ + ce-^+W] 



(39) 



where b and c are complex constants determined by the initial conditions. Stable modes 
will be described by (3 taking a real nonintegral value (we exclude the case of integral /3). 
Following the discussion in §4.70 of 66|, for trapping parameters in the first stability zone 
of the Mathieu equation, /3 can be chosen in the range < (3 < 1 for all stable modes. 
By defining i?2n = A — (2n + f3) 2 we can write infinite recursion relations for C2 n , 

QC2n-2 — R2nC2n ~ QC2n+2i (40) 

and obtain two independent expansions in infinite continued matrix inversions 



C 2 = T 2 ,pQC = ( \R 2 - Q [R 4 - Q [R 6 - ...p 1 Q] 1 Q 



-1 



QC . 



and 



QC 2 = RqCq — QC-2 = Tq^Cq = (r — Q [R-2 — Q 4 — •••] 1 Q] 0] Co 



Multiplying eq. f )4T|) by Q and defining 



^2,/3 = Tq$ — QT 2 ^Q, 



(41) 



(42) 



(43) 



we find that all characteristic exponents (3 are zeros of the determinant of Y 2j p (which is 
a function of 0). If there are degenerate /3's they will appear as degenerate zeros of this 
determinant. The vector Co for each (3 is an eigenvector of Y 2 ^ with eigenvalue 0. Since A 
and Q are symmetric, Y 2y p is symmetric as well, and so its kernel will be of dimension equal 
to the algebraic multiplicity of the (3 root. The vector C2 can be obtained by an application 
of T 2 ^Q to Co, for n = — 1 we use (7_ 2 = [T^,^ 1 QCo, an d so on for the other vectors. We 
note that the different vectors Cm$ are not orthogonal in general, and the vectors at every 
order in n mix different coordinates. 

The general term of the expansion vanishes. Either A or Q may be singular and the 
expansion can still be applied in general. Even if both are singular, the expansion is valid if 
there are no integral values of /3, a case which we do not tackle as noted above. Excluding 
perhaps isolated values of (3 (and atypically in the a — q parameter-space), all matrices which 
are inverted in the above expressions will be invertible, and while employing the algorithm 
in practice, the invertibility of the matrices is of course easily verified at each step. In Sec. [V] 
we use in fact a generalization of the above expansion [72| which includes also the next 
Fourier harmonic (cos4t, ommited from eq. fl30|) ). 



C. The Floquet-Lyapunov Transformation For Stable Modes 



We now further assume that all Floquet modes are stable, i.e. that the 2f linearly- 
independent solutions of eq. (J3"4"}) are oscillatory and thus come in complex conjugate pairs. 
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This simplifies many expressions. We therefore take B of eq. ( |35|) in the block form 

B = -is)' (^)/x/ = diag{/3 1 ,.. (44) 

where /3j are positive. We define the /-dimensional matrix U whose columns are constructed 
from the series of /-dimensional vectors C^p. obtained from the recursion relations for the 
solutions of eq. (|39|) . i.e. 

(P)fxf= (E^2n,^ 2nt -) (45) 

where in the above expression and for the rest of this section, the summation is over iigZ. 
We similarly define the /-dimensional matrix V composed from column- vectors as 

QOfxf = ( i E (2n + Pj) C 2n ^ nt ... ) . (46) 
The matrices U and V can be chosen to obey the normalization condition 

V* (0) U (0) = \i, (47) 
which is a rescaling imposed by multiplication with a (diagonal) matrix, such that 

U->U {-21V 1 (0) C/ (0)) ~* , (48) 



and V accordingly. As shown in 72|], the Floquet-Lyapunov transformation and its inverse 
can then be obtained in closed form, and are given block-wise by (where U* denotes the 
complex conjugate of the matrix U), 

This transformation is a canonical transformation from the Hamiltonian coordinates u 
and their conjugate momenta p — u, to the variables given by 



X 



where £ is the new /-dimensional vector of coordinates and — i( are the new conjugate 
momenta (we here break the notation a little). Using the realness of u and p it is easy to 
verify that £ = (*. The time dependence of these modes is 



£i (t) = $ (0) 0(t)=0(0)e~^. (50) 



V. COMPUTATION OF THE MODES OF A 'PECULIAR' CRYSTAL 



In Figs. [T]|3]we consider a crystal of 6 ions, demonstrating the utility and generality of the 
analysis presented in this contribution. The simulation is of an almost ideal linear Paul trap 



15 




N 



N 



FIG. 1. The periodic 'peculiar' crystal solution of 6 ions in a linear Paul at nearly-spherical 
trapping parameters (see text for details). The average locations correspond to a (nearly-regular) 
Octahedron. Upper-left : the ions at zero rf-phase. The other 3 figures show the ion trajectories 
over one rf-period. The oscillation with a large amplitude is given by eq. (|23p with an accuracy of 
0.5%. The axial micromotion amplitude is of order 10 -4 of the ions' radial positions (see Sec. lIII Cj) . 



(with only a 1% DC asymmetry in the radial plane), such that the center-of-mass frequencies 
are degenerate to within 1%. The Mathieu parameters are q y = 0.41 and a x = 0.05766. In 
the pseudopotential approximation, the corresponding minimum-configuration is a (nearly- 
regular) octahedron, with two ions sitting on each axis of the trap. The rf crystal on the 
other hand, may well deserve to be called 'peculiar' (in the sense described in Sec. [TTJ) , since 
it is not oriented with the axes. There are two ions lying along the y-axis, along which the 
confinement is strongest, but the two other pairs of ions sit along lines which are rotated in 
x — z plane. The large amplitude oscillation at the rf-period, and the absence of micromotion 
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FIG. 2. Ions' primary direction of oscillation (given by Co), in four Floquet-Lyapunov modes, 
which are numbered from the lowest frequency, with £1 being the lowest-frequency mode. 



along the axial direction, confirm very accurately the results presented in Sec. IIII CI 

Numerically, the periodic crystal solution can be obtained by starting from a simulation 
of the full e.o.m's with a friction (cooling) term and slowly turning it off (the adiabatic 
shutting-down of the damping term is important). The crystal is then followed for a period 
and the periodic solution and force matrix can be Fourier expanded to obtain the matrices A 
and Q. Since the crystal is 'peculiar', it has no corresponding pseudopotential limit (normal 



respc 

modes of regular polyhderons were investigated in 751 ]. and see also references in [76j). For a 
very accurate description of the modes in the Paul trap, we add to eq. ( 132]) the next Fourier 
harmonic to get u + [A — 2Q 2 cos2i — 2Q 4 cos4t] u — 0, and expand it in the method of 
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FIG. 3. Comparison between the analytical expansion and the exact numerical solution, for the 
Floquet-Lyapunov modes of Fig. [21 with small random initial conditions. The modes are numbered 
from the lowest frequency, with £i being the lowest-frequency mode. Time is measured in periods 
of the rf frequency, and corresponding natural nondimensional units are used for distances. 



continued matrix inversions, using the formulas given in App. B of 72j . which generalize the 
formulas presented in Sec. IIVI This modification is required in order to obtain accurately 
the low frequency modes of nearly-degenerate configurations, as in a peculiar crystal. 
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VI. CONCLUDING COMMENTS 



In this contribution we have investigated the dynamics of ion crystals in rf traps. We 
repeat the main results presented herein. In eq. ([7]) we show that in general, the rf couples 
the pseudopotential normal modes of the crystal, and also acts as a driving force. For a 
crystal configuration which has the same symmetry as the trapping potential (e.g. a chain 
of ions or a planar crystal in a trap with cylindrical symmetry), we find that the equations 
of motion for the modes become decoupled Mathieu equations as in eq. f JTOj) . for which the 
stability analysis is trivial (but may be different from the pseudopotential approximation). 

In eqs. fl23|) and f l27|) we derive results which have been observed in experiments and 
numerical simulations, namely that in a general ion crystal the micromotion amplitude in 
each coordinate is q a /2 of the respective average position, and that in the linear Paul trap, 
axial micromotion is negligibly small. 

When the crystal solution differs from the pseudopotential limit (as in the peculiar crystals 
discussed above), or when a very accurate expansion of the modes is desired, the derivation of 
eq. ( 1321) takes into account the full rf crystal solution and ion interactions along the periodic 
trajectory. Sec. HV] describes briefly how to solve for the decoupled modes of oscillations of 
the ions, allowing obtaining explicitly the mode frequencies and solution vectors. 

We have focused on single-species crystals in quadrupole traps, specifically the linear and 
hyperbolic Paul traps. However, generalization to other cases is easy. Crystals in multipole 
traps 77-7^| can be treated by expanding the motion around a suitable periodic solution 
RJ \ as in Sec. II II CI Keeping only the leading terms will lead to the equations treated 



above. Segmented traps and trap arrays [80h83| can be handled simply by changing the 
Mathieu parameters felt by each ion (assuming that the rf drive has identical frequency). 
Crystals of nonidentical ions and other types of driving can be treated similarly, and various 
transformations [84| can be used to handle more general linear system similar to eq. (|32|) . such 
as systems with first-order derivatives (e.g. linear damping, gyroscopic forces or magnetic 
fields BH)- 

The framework presented here for calculating the classical normal modes and their fre- 
quencies can find immediate application in most studies involving trapped ion Coulomb 
crystals .^including studies of how ion Coulomb crystals differ from uniformly charged liquid 
models [I4j, and possibly even crystalline beams 86|, |87J. The analysis is based on lineariza- 
tion of the nonlinear solution about a periodic solution, and the nonlinear correction terms 
can be written as a series expansion in the Floquet modes about the crystal solution. In a 
quadrupole trap the nonlinear terms would come only from the Coulomb interaction, and 
in a multipole trap, there will be terms coming from the nonlinear trapping potential. Such 
an expansion may serve as a starting point for study of nonlinear collective phenomena, e.g. 



20, |73|, |8£l |89| and many more) 



the important phenomenon of rf heating (e.g. 

In addition, an exact quantum description of the modes, and wavefunctions in the configu- 
ration space of the ions, are presented in 72|, utilizing the Floquet-Lyapunov transformation 
described in Sec. [IVl The nontrivial time-dependent wavefunctions could, e.g., eventually 
become important for understanding trapped ion chemistry at ultracold temperatures 90 . 
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